#include "mex.h"
#include "matrix.h"
#include <cstdio>
#include <cstdlib>
#include <vector> 
#include <time.h>
#include <errno.h>
#include "./source/pnmfile.h"
#include "./source/filter.h"
#include "./source/image.h"
#include "./source/misc.h"
#include "./source/matrix.h"

    //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
    //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
    //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

using namespace vlib;
using namespace std;

double large = 1E10F;

// error checking
#define check(e) \
(e ? (void)0 : (printf("%s:%u error: %s\n%s\n", __FILE__, __LINE__, #e, strerror(errno)), exit(1)))

// compute integral along each column of an image
// I(x, y) = integral of column x, from row 0 to row y-1
// the -1 offset helps simplyfing the DP code
image<double> *col_integral(image<double> *im) {
  int width = im->width();
  int height = im->height();
  image<double> *I = new image<double>(width, height+1);
  for (int x = 0; x < width; x++) {
    imRef(I, x, 0) = 0;
    for (int y = 1; y <= height; y++) {
      imRef(I, x, y) = imRef(I, x, y-1) + imRef(im, x, y-1);
    }
  }
  return I;
}

// These macros are used to access the DP tables
#define Pref(k,l,i,j)						\
  (((((k) * K) + (l)) * (height+1) + (i)) * (height+1) + (j))
#define Vref(l,i,j)				\
  (((l) * (height+1) + (i)) * (height+1) + (j))
// These macros are used pack a state into a single integer
#define pack(l,i,j)				\
  ((l) | ((i) << 10) | ((j) << 20))
#define unpack(P)	       \
  int l = (P) & 0x3FF;         \
  int i = ((P) >> 10) & 0x3FF; \
  int j = ((P) >> 20) & 0x3FF;			

// init for propagation
#define initprop()							\
  for (int i = 0; i <= height; i++) {					\
    for (int j = i; j <= height; j++) {					\
      double disc = imRef(WH, k-1, i) * matRef(horz, 0,0) +		\
	(imRef(WH, k-1, height)-imRef(WH, k-1, j)) * matRef(horz,1,1) +	\
	(imRef(WH, k-1, j)-imRef(WH, k-1, i)) * matRef(horz,a,b);	\
      Vtmp[Vref(0,i,j)] = Vold[Vref(a,i,j)] + disc;			\
      Ptmp[Vref(0,i,j)] = pack(a,i,j);					\
    }									\
  }

// propagate increasing i to j
#define increasei()							\
  for (int j = 1; j <= height; j++) {					\
    for (int i = 1; i < j; i++) {					\
      double v = Vtmp[Vref(0,i-1,j)] + imRef(whorz, k-1, i-1) * (matRef(horz,a,0) - matRef(horz,a,b)); \
      if (v < Vtmp[Vref(0,i,j)]) {					\
	Vtmp[Vref(0,i,j)] = v;						\
	Ptmp[Vref(0,i,j)] = Ptmp[Vref(0,i-1,j)];			\
      }									\
    }									\
  }

// propagate decreasing j to i
#define decreasej()							\
  for (int i = 0; i < height; i++) {					\
    for (int j = height-1; j > i; j--) {				\
      double v = Vtmp[Vref(0,i,j+1)] + imRef(whorz, k-1, j) * (matRef(horz,a,1) - matRef(horz,a,b)); \
      if (v < Vtmp[Vref(0,i,j)]) {					\
	Vtmp[Vref(0,i,j)] = v;						\
	Ptmp[Vref(0,i,j)] = Ptmp[Vref(0,i,j+1)];			\
      }									\
    }									\
  }								

// propagate decreasing i
#define decreasei()						\
  for (int j = 1; j <= height; j++) {				\
    for (int i = j-1; i >= 0; i--) {				\
      double v = Vtmp[Vref(0,i+1,j)] + imRef(whorz, k-1, i) * (matRef(horz,0,b) - matRef(horz,0,0)); \
      if (v < Vtmp[Vref(0,i,j)]) {				\
	Vtmp[Vref(0,i,j)] = v;					\
	Ptmp[Vref(0,i,j)] = Ptmp[Vref(0,i+1,j)];		\
      }								\
    }								\
  }

// propagate increasing j
#define increasej()						\
  for (int i = 0; i < height; i++) {				\
    for (int j = i+1; j <= height; j++) {			\
      double v = Vtmp[Vref(0,i,j-1)] + imRef(whorz, k-1, j-1) * (matRef(horz,1,b) - matRef(horz,1,1)); \
      if (v < Vtmp[Vref(0,i,j)]) {				\
	Vtmp[Vref(0,i,j)] = v;					\
	Ptmp[Vref(0,i,j)] = Ptmp[Vref(0,i,j-1)];		\
      }								\
    }								\
  }	

// update table 
#define updateprop()				      \
  for (int i = 0; i <= height; i++) {		      \
    for (int j = i; j <= height; j++) {		      \
      if (Vtmp[Vref(0,i,j)] < V[Vref(b,i,j)]) {	      \
	V[Vref(b,i,j)] = Vtmp[Vref(0,i,j)];	      \
	P[Pref(k,b,i,j)] = Ptmp[Vref(0,i,j)];	      \
      }						      \
    }						      \
  }


/*
 * Labeling algorithm
 *
 * image<uchar> *label(int K, image<double> **data, 
 *                     matrix<double> *horz, matrix<double> *vert,
 *		       image<double> *whorz, image<double> *wvert)
 *
 * The function returns an optimal tiered labeling
 *
 * k is the number of labels, must be at least 3
 * 0 is the top region label
 * 1 is the bottom region label
 * 2 ... (k-1) are the middle region labels
 *
 * data is an array of k images, defining data costs for each label
 *
 * horz, vert are matricies with the discontinuity costs
 * matRef(horz, a, b) is the cost of having a on the left of b
 * matRef(vert, a, b) is the cost of having a on top of b
 * 
 * whorz, wvert are spatially varying discontinuity costs
 * imRef(whorz, x, y) is the cost between (x,y) and (x+1,y)
 * imRef(wvert, x, y) is the cost between (x,y) and (x,y+1)
 *
 * The two costs are multiplied together
 * If pixel (x,y) is labeled a, and pixel (x+1,y) is labeled b the 
 * discontinuity cost is:
 * matRef(horz, a, b) * imRef(whorz, x, y)
 *
 * NOTE: This version enforces the middle region in column k to intersect 
 * the middle region in column k-1, unless one of those regions is empty. 
 * This is typically what we want, and eliminates 2 cases.
 */
image<uchar> *label(int K, image<double> **data, 
		    matrix<double> *horz, matrix<double> *vert,
		    image<double> *whorz, image<double> *wvert) {
  int width = data[0]->width();
  int height = data[0]->height();
  check(width < (1 << 11));
  check(height < (1 << 11));
  check(K < (1 << 11));

  // column integrals of data costs
  vector<image<double> *> I(K);

  for (int i = 0; i < K; i++)
    I[i] = col_integral(data[i]);

  image<double > *WH, *WV;
  WH = col_integral(whorz);
  WV = col_integral(wvert);
  
  /* 
   * DP table 
   * 
   * V[l,i,j] is the cost of an optimal labeling up to (and
   * including) the current column such that the top region ends in row (i-1),
   * the bottom region starts in row j and l is the label of the
   * middle region.
   *
   * The end of the top region is defined to be (i-1) so that when i=0
   * the top region is empty.  If j = i there is no middle region.
   * Entries with j < i are undefined.
   *
   * We keep traceback pointers in P[k,l,i,j] for each column k.
   */

  double *V = new double[K*(height+1)*(height+1)];
  int *P = new int[width*K*(height+1)*(height+1)];
  double *Vold = new double[K*(height+1)*(height+1)];
  double *Vtmp = new double[(height+1)*(height+1)];
  int *Ptmp = new int[(height+1)*(height+1)];

  memset(V, 0, K*(height+1)*(height+1)*sizeof(double));
  for (int k = 0; k < width; k++) {
    memcpy(Vold, V, K*(height+1)*(height+1)*sizeof(double));
    
    // for each label of the middle region in this column
    for (int b = 2; b < K; b++) {
      if (k > 0) {
	for (int i = 0; i <= height; i++) {
	  for (int j = i; j <= height; j++) {
	    V[Vref(b,i,j)] = INF;
	  }
	}
	
	// for each label of the middle region in the previous column
	for (int a = 2; a < K; a++) {
	  initprop();
	  increasei();
	  decreasej();
	  updateprop();
	  
	  initprop();
	  increasei();
	  increasej();
	  updateprop();
	  
	  initprop();
	  decreasej();
	  decreasei();
	  updateprop();
	  
	  initprop();
	  increasej();
	  decreasei();
	  updateprop();
	}
	
	// allow the top/bottom boundary to move when there is no middle region
	for (int i = 0; i <= height; i++) {
	  for (int ip = 0; ip <= height; ip++) {
	    double h = 0;
	    if (ip < i)
	      h = imRef(WH, k-1, ip) * matRef(horz, 0, 0) +
		(imRef(WH, k-1, i)-imRef(WH, k-1, ip)) * matRef(horz, 1, 0) +
		(imRef(WH, k-1, height)-imRef(WH, k-1, i)) * matRef(horz,1,1);
	    else 
	      h = imRef(WH, k-1, i) * matRef(horz, 0, 0) +
		(imRef(WH, k-1, ip)-imRef(WH, k-1, i)) * matRef(horz, 1, 0) +
		(imRef(WH, k-1, height)-imRef(WH, k-1, ip)) * matRef(horz,1,1);
	    double v = Vold[Vref(b,ip,ip)] + h;
	    if (v < V[Vref(b,i,i)]) {
	      V[Vref(b,i,i)] = v;
	      P[Pref(k,b,i,i)] = pack(b,ip,ip);
	    }
	  }
	}
      }

      // add in data term and vertical discontinuities
      for (int i = 0; i <= height; i++) {
	for (int j = i; j <= height; j++) {
	  double data_top = imRef(I[0], k, i);
	  double data_bottom = imRef(I[1], k, height) - imRef(I[1], k, j);
	  double data_middle = imRef(I[b], k, j) - imRef(I[b], k, i);
	  double disc = 0;
	  // top middle boundary
	  if (i > 0 && j > i)
	    disc += matRef(vert, 0, b) * imRef(wvert, k, i-1);
	  // middle bottom boundary
	  if (j > i && j < height)
	    disc += matRef(vert, b, 1) * imRef(wvert, k, j-1);
	  // top bottom boundary
	  if (i > 0 && j < height && i == j)
	    disc += matRef(vert, 0, 1) * imRef(wvert, k, i-1);
	  // vertical discontinuities within top
	  if (i > 1)
	    disc += matRef(vert, 0, 0) * imRef(WV, k, i-1);
	  // vertical discontinuities within bottom
	  if (j < height-1)
	    disc += matRef(vert, 1, 1) * (imRef(WV, k, height) - imRef(WV, k, j));
	  // vertical discontinuities within middle
	  if (j > i)
	    disc += matRef(vert, b, b) * (imRef(WV, k, j-1) - imRef(WV, k, i));
	  V[Vref(b,i,j)] += data_top + data_bottom + data_middle + disc;
	}
      }
    }
  }

  // trace back
  image<uchar> *out = new image<uchar>(width, height);  

  double best = INF;
  int bi, bj, bl;
  for (int l = 2; l < K; l++) { 
    for (int i = 0; i <= height; i++) {
      for (int j = i; j <= height; j++) {
	double val = V[Vref(l, i, j)];
	if (val < best) {
	  best = val;
	  bi = i;
	  bj = j;
	  bl = l;
	}
      }
    }
  }
  //printf("best = %f\n", best);
  int x = width-1;
  while (x >= 0) {
    for (int y = 0; y < height; y++) {
      if (y < bi)
	imRef(out, x, y) = 0;
      else if (y >= bj)
	imRef(out, x, y) = 1;
      else
	imRef(out, x, y) = bl;
    }
    int p = P[Pref(x, bl, bi, bj)];
    unpack(p);
    bl = l;
    bi = i;
    bj = j;
    x--;
  }

  for (int i = 0; i < K; i++)
    delete I[i];
  delete WH;
  delete WV;
  delete V;
  delete P;
  delete Vold;
  delete Vtmp;
  delete Ptmp;

  return out;
}

/*
 * Compute energy of labeling 
 */
double energy(image<uchar> *l, image<double> **data,
	      matrix<double> *horz, matrix<double> *vert, 
	      image<double> *whorz, image<double> *wvert) {
  int width = l->width();
  int height = l->height();
  
  double e = 0;
  for (int y = 0; y < height; y++) {
    for (int x = 0; x < width; x++) {
      int a = imRef(l, x, y);
      e += imRef(data[a], x, y);
      if (x < width-1) {
	int b = imRef(l, x+1, y);
	e += matRef(horz, a, b) * imRef(whorz, x, y);
      }
      if (y < height-1) {
	int b = imRef(l, x, y+1);
	e += matRef(vert, a, b) * imRef(wvert, x, y);
	
      }
    }
  }
  return e;
}

/*
 * Find optimal labeling with a single connected middle/foreground region
 *
 * This requires the top and bottom labels use the same data costs
 *
 * We reduce the problem to the standard tiered labeling problem 
 * The trick: there are two fresh "dummy" middle labels, one indicates the 
 * foreground region hasn't started,
 * the other indicates the foreground region has ended.  
 *
 */
image<uchar> *label_con(int K, image<double> **data,
			matrix<double> *horz, matrix<double> *vert,
			image<double> *whorz, image<double> *wvert) {
  // data for K and K+1 is same as top
  image<double> **data_con = new image<double> *[K+2];
  for (int i = 0; i < K; i++)
    data_con[i] = data[i];
  data_con[K] = data[0];
  data_con[K+1] = data[0];

  matrix<double> *horz_con = new matrix<double>(K+2, K+2);
  matrix<double> *vert_con = new matrix<double>(K+2, K+2);
  horz_con->init(large);
  vert_con->init(large);
  for (int i = 0; i < K; i++) {
    for (int j = 0; j < K; j++) {
      matRef(horz_con, i, j) = matRef(horz, i, j);
      matRef(vert_con, i, j) = matRef(vert, i, j);
    }
  }

  // make top/bottom boundary forbidden
  // instead we must stick K or K+1 in between
  matRef(vert_con, 0, 1) = large; 
  matRef(vert_con, 0, K) = matRef(vert, 0, 0);
  matRef(vert_con, 0, K+1) = matRef(vert, 0, 0);
  matRef(vert_con, K, 1) = matRef(vert, 0, 1);
  matRef(vert_con, K+1, 1) = matRef(vert, 0, 1);
  
  // we use K until the foreground region, afterwards we use K+1
  matRef(horz_con, K, K) = 0;
  matRef(horz_con, K+1, K+1) = 0;
  for (int i = 2; i < K; i++) {
    matRef(horz_con, K, i) = matRef(horz, 0, i);
    matRef(horz_con, K+1, i) = large;
    matRef(horz_con, i, K) = large;
    matRef(horz_con, i, K+1) = matRef(horz, i, 0);
  }
  
  matRef(horz_con, K, 0) = large;
  matRef(horz_con, K, 1) = large;
  matRef(horz_con, 0, K) = large;
  matRef(horz_con, 1, K) = large;
  matRef(horz_con, K+1, 0) = large;
  matRef(horz_con, K+1, 1) = large;
  matRef(horz_con, 0, K+1) = large;
  matRef(horz_con, 1, K+1) = large;
  
  image<uchar> *out = label(K+2, data_con, horz_con, vert_con, whorz, wvert);
  
  int width = out->width();
  int height = out->height();
  for (int y = 0; y < height; y++) {
    for (int x = 0; x < width; x++) {
      if (imRef(out, x, y) >= K)
	imRef(out, x, y) = 0;
    }
  }

  delete horz_con;
  delete vert_con;
  delete data_con;
  return out;
}

/*
 * Save segmentation illustrating part decomposition
 */
void saveSegmentation(image<uchar> *segmentation, int K, const char *name) {
  int scale = 255/K;
  int width = segmentation->width();
  int height = segmentation->height();
  image<uchar> *tmp = new image<uchar>(width, height);

  for (int y = 0; y < height; y++) {   
    for (int x = 0; x < width; x++) {
      /*if (imRef(segmentation,x,y) < 2)
	imRef(tmp,x,y) = 255;
      else*/ imRef(tmp,x,y) = imRef(segmentation,x,y)*scale;
    }
  }

  savePGM(tmp, name);
  delete tmp;
}

/*
 * Simple data model for foreground/background segmentation
 * We assume the foreground is bright and background is dark
 */
void loadData(image<uchar> *I, image<double> * &fg, image<double> * &bg) {
  int width = I->width();
  int height = I->height();
  
  fg = new image<double>(width, height);
  bg = new image<double>(width, height);
  double robustThreshold = 300;//30; dolev
  for (int y = 0; y < height; y++) {
    for (int x = 0; x < width; x++) {
      imRef(fg, x, y) = abs(150-imRef(I, x, y));
      imRef(bg, x, y) = abs(50-imRef(I, x, y));
      if (imRef(fg, x, y ) > robustThreshold) 
	imRef(fg, x, y ) = robustThreshold;
      if (imRef(bg, x, y ) > robustThreshold) 
	imRef(bg, x, y ) = robustThreshold;
    }
  }
}

/*
 * Simple discontinuity costs based on local contrast
 */
void initializeW(image<uchar> *I, 
		 image<double> * &whorz,image<double> * &wvert) {  
  int width = I->width();
  int height = I->height();

  double aver = 0.0;
  int total = 0;
  for (int y = 0; y < height-1; y++) {
    for (int  x = 0; x < width-1; x++) {
      aver = aver + 
	abs(imRef(I,x,y)-imRef(I,x+1,y))+
	abs(imRef(I,x,y)-imRef(I,x,y+1));
      total = total+2;
    }
  }
  aver = aver/total;	

  whorz = new image<double>(width, height);
  wvert = new image<double>(width, height);

  for (int y = 0; y < height; y++) {
    for (int x = 0; x < width-1; x++) {
      double diff = abs(imRef(I,x+1,y)-imRef(I,x,y));
      if (diff > 3.0*aver) 
	imRef(whorz,x,y) = 1;
      else if (diff > 2.0*aver) 
	imRef(whorz,x,y) = 2;
      else if (diff > 1.5*aver)
	imRef(whorz,x,y) = 3;
      else imRef(whorz,x,y) = 8;
    }
  }

  for (int  y = 0; y < height-1; y++) {
    for (int x = 0; x < width; x++ ) {
      double diff = abs(imRef(I,x,y+1)-imRef(I,x,y));
      if (diff > 3.0*aver) 
	imRef(wvert,x,y) = 1;
      else if (diff > 2.0*aver) 
	imRef(wvert,x,y) = 2;
      else if (diff > 1.5*aver) 
	imRef(wvert,x,y) = 3;
      else imRef(wvert,x,y) = 8;
    }
  }
}
 

/*
 * Initialize data and shape model
 */
void loadModel(image<uchar> *I, int &K, 
	       matrix<double> * &horz, matrix<double> * &vert,
	       image<double> ** &data) {
  K = 4;
  horz = new matrix<double>(K,K);
  vert = new matrix<double>(K,K);
    
  double lambda = 20;
  horz->init(lambda);
  vert->init(lambda);

  for (int i = 0; i < K; i++ ) {
    matRef(horz,i,i) = 0;
    matRef(vert,i,i) = 0;
  }
  
  int inMatHorz[4][4] = {	  {0,1000,5,100},
						  {1000,0,100,100},
						  {1000,100,0,5},
						  {5,100,1000,0}       };

  int inMatVert[4][4] = {	  {0,100,5,5},
						  {1000,0,1000,1000},
						  {1000,5,0,100},
						  {1000,5,100,0}       };

  for (int i_index = 0; i_index < K ; i_index++) {
	for (int j_index = 0; j_index < K ; j_index++) {
	  matRef(horz,i_index,j_index) = inMatHorz[i_index][j_index] * lambda;
	  matRef(vert,i_index,j_index) = inMatVert[i_index][j_index] * lambda;
	}
  }
  /*
  int num = 4;
  // top - 0
  // midle - 1
  // bottom - 2
  matRef(horz,0,1) = 160 * lambda;
  matRef(horz,1,0) = num * lambda;

  matRef(horz,2,1) = num * lambda;
  matRef(horz,1,2) = num * lambda;

  matRef(horz,0,2) = large;
  matRef(horz,2,0) = large;

  num = 1;
  matRef(vert,0,1) = num * lambda;
  matRef(vert,0,2) = large;

  matRef(vert,1,0) = num * lambda;
  matRef(vert,1,2) = large;

  matRef(vert,2,0) = large;
  matRef(vert,2,1) = large;
*/


/*  for (int i = 0; i < K; i++ ) {
    matRef(horz,i,i) = 0;
    matRef(vert,i,i) = 0;
  }

  // 0/1 are both background
  // so there is no discontinuity cost between them
  matRef(vert,0,1) = lambda;
  matRef(vert,1,0) = large;
  matRef(horz,0,1) = lambda*5;
  matRef(horz,1,0) = lambda*5;
	
  // 2/3 are both foreground
  // so there is no discontinuity cost between them, 
  // except that 2 is the left side, 3 is the right side
  matRef(horz,2,3) = 0;
  matRef(horz,3,2) = 0;  // R,L transition forbidden	

  // we forbid transitions to enforce a slant of the boundary
//  matRef(horz,2,0) = 5*lamda;  // L,T transition
//  matRef(horz,2,1) = 5*lamda;  // L,B transition
  matRef(horz,0,3) = large;  // T,R transition
  matRef(horz,1,3) = large;  // B,R transition
	
  // we discourage vertical transitions more than horizontal
  // so that the boundaries are more vertical
  int mult = 1;
  matRef(vert,0,2) = lambda;
  matRef(vert,2,1) = lambda;
  matRef(vert,0,3) = lambda;
  matRef(vert,3,1) = lambda;  
  */
  // data model
  image<double> *fg, *bg;
  loadData(I, fg, bg);
  data = new image<double> *[K];
  data[0] = bg;
  data[1] = bg;
  for (int i = 2; i < K; i++)
    data[i] = fg;
}
    //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
    //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
    //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
    
void finishSegmantation(image<uchar> *segmentation, int K) {
      int scale = 255/K;
      int width = segmentation->width();
      int height = segmentation->height();

      for (int y = 0; y < height; y++) {   
        for (int x = 0; x < width; x++) {
           imRef(segmentation,x,y) *= scale;
        }
      }
}

image<uchar>* maintired(image<uchar> *input) {
    //this method is a replication of the main function in tiered labeling

    /////////////////////////////////////////////////////////////////////////////////////////////
    // copied code from tiered labeling main function    //
    /////////////////////////////////////////////////////////////////////////////////////////////
      int K;
      vlib::matrix<double> *horz,  *vert;
      vlib::image<double>  *whorz, *wvert;
      vlib::image<double>  **data;
      
      // the load pgm
      
      // initialize spatially varying discontinuity costs
     initializeW(input, whorz, wvert);

      // initialize shape and data model
      loadModel(input, K, horz, vert, data);

      // find optimal labeling
      // forcing middle region to be one connected component
    image<uchar> *l = label(K, data, horz, vert, whorz, wvert);
    /////////////////////////////////////////////////////////////////////////////////////////////
    // finished the copied section                                                  //
    /////////////////////////////////////////////////////////////////////////////////////////////
    
    // overiden method of save segmantaion
    finishSegmantation(l,  K);
    
    delete input;
    delete horz;
    delete vert;
    delete whorz;
    delete wvert;
    delete data;
    
    return l;
    return NULL;
}


void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
    // test args
    /* check for proper number of arguments */
    if(nrhs!=1) {
        mexErrMsgIdAndTxt("MyToolbox:arrayProduct:nrhs", "One inputs required.");
    }
    if(nlhs!=1) {
        mexErrMsgIdAndTxt("MyToolbox:arrayProduct:nlhs",  "One output required.");
    }
    
    /* make sure the input argument is image */
    if( !mxIsDouble(prhs[0])) {
        mexErrMsgIdAndTxt("MyToolbox:arrayProduct:prhs", "Input must be uint8 matrix.");
    }
    
    // convert input to image
    int width = mxGetM (prhs[0]);
    int height = mxGetN (prhs[0]);
    vlib::image<vlib::uchar> *input = new vlib::image<vlib::uchar>(width, height);
    
    // insert data to image objec (equivalent to loadPGM()  )
    double *inMatrix;
    inMatrix = mxGetPr(prhs[0]);
    for (int x=0; x < width; x++) {
       for (int y=0; y < height; y++) {
           imRef(input, y, x) = inMatrix[y*width+x];
        }
    }
    
    // calling the tired labling logic
    vlib::image<vlib::uchar> *output = maintired(input);
    
  
    /* create the output matrix */
    plhs[0] = mxCreateDoubleMatrix(height,width,mxREAL);

    /* get a pointer to the real data in the output matrix */
    double *outMatrix = mxGetPr(plhs[0]);
    // populate matrix
    for (int x=0; x < width; x++) {
        for (int y=0; y < height; y++) {
            outMatrix[y*width+x] = imRef(output, y, x);
        }
    }
    delete output;
}